% !TEX root = FDS_Technical_Reference_Guide.tex

\typeout{new file: Mass_Chapter.tex}

\chapter{Mass, Species, and Enthalpy Transport}

This chapter describes in detail the equation of state in the low Mach number limit, the finite difference approximation of the mass and species conservation equations, and the role of the flow divergence as a surrogate for the enthalpy transport equation.
Due to the use of the low Mach number approximation, the energy conservation equation is not solved explicitly but rather is defined implicitly via the divergence of the flow field, which contains the combustion and radiation source terms.


\section{The Equation of State}

A distinguishing feature of a CFD model is the regime of
flow speeds (relative to the speed of sound) for which it is designed. High
speed flow codes involve compressibility effects and shock waves. Low speed
solvers, however, explicitly eliminate compressibility effects that give rise
to acoustic (sound) waves. The Navier-Stokes equations describe the
propagation of information at speeds comparable to that of the fluid flow (for fire, approximately \SI{10}{m/s}),
but also at speeds comparable to that of sound waves (for still air,
\SI{300}{m/s}). Solving a discretized form of these equations would require extremely small
time steps in order to account for information traveling at the speed of sound, making
practical simulations difficult.

Following the work of Rehm and Baum~\cite{Rehm:1}, an approximation to the equation of state is made by decomposing the pressure
into a ``background'' component and a perturbation. It is assumed that
the background component of pressure can differ from compartment to compartment. If
a volume within the computational domain is isolated from other volumes, except via leak paths or ventilation ducts, it is referred to as a ``pressure zone'' and assigned its own background pressure. The pressure field within the $m$th zone, for example, is a linear combination
of its background component and the flow-induced perturbation:
\be
   p(\bx,t) = \bp_m(z,t) + \tp(\bx,t)
\ee
Note that the background pressure is a function of $z$, the vertical spatial coordinate, and the time, $t$. For most compartment fire applications, $\bp_m$ changes very little with height or time. However, for scenarios where a fire increases the pressure in a closed compartment, or where the HVAC system affects the pressure, or when the height of the domain is significant, $\bp_m$ takes these effects into account~\cite{Baum:5}. The ambient pressure field is denoted $\bp_0(z)$. Note that the subscript 0 denotes the exterior of the computational domain, not time 0. This is the assumed atmospheric pressure stratification that serves as both the initial and boundary condition for the governing equations.

The purpose of decomposing the pressure is that for low Mach number flows, it can be assumed that the temperature and density are inversely proportional, and thus the equation of state (in the $m$th pressure zone) can be approximated as
\be
   \bp_m  =  \rho T \R \sum_\alpha \frac{Z_\alpha}{W_\alpha} = \frac{\rho T \R}{ \bW }  \label{state}
\ee
Recall from Section~\ref{sec_lumped_species} that $Z_\alpha$ is the mass fraction of lumped species $\alpha$. The pressure, $p$, in the state and energy equations is replaced by the background pressure $\bp_m$ to filter out sound waves that travel at speeds that are much faster than typical flow speeds expected in fire applications. The low Mach number assumption serves two purposes. First, the filtering of acoustic waves means that the time step in the numerical algorithm is bound only by the flow speed as opposed to the speed of sound, and second, the modified state equation leads to a reduction in the number of dependent variables in the system of equations by one. The energy equation (\ref{energy}) is not explicitly solved; rather, its source terms are included in the expression for the flow divergence, to be discussed later in the chapter.  When the velocity field satisfies the specified thermodynamic divergence, the conservative form of the sensible enthalpy equation is satisfied by construction.

The stratification of the atmosphere is derived from the relation
\be \frac{\d \bp_0}{\d z} = - \rho_0(z) \, g  \ee
where $\rho_0$ is the background density and $g=\SI{9.8}{m/s}$. Using Eq.~(\ref{state}), the background pressure can be written as a function of the background temperature, $T_0(z)$,
\be \bp_0(z) = p_\infty \; \exp \, \left( -\int^z_{z_\infty} \frac{\bW \, g}{\R \, T_0(z')} \d z' \right)  \label{pstrat} \ee
where the subscript infinity generally refers to the ground. A linear temperature stratification of the atmosphere may be
specified by the user such that $T_0(z) = T_\infty + \Gamma z$ where $T_\infty$ is the temperature at the ground and
$\Gamma$ is the lapse rate (e.g., $\Gamma = -\SI{0.0098}{K/m}$ is the {\em adiabatic lapse rate}).
In this case $\bp_0$ and $\rho_0$ are derived from Eqs.~(\ref{pstrat}) and (\ref{state}), respectively.
It can then be shown that for $\Gamma \ne 0$ the pressure stratification becomes
\be
   \bp_0(z) = p_\infty  \left( \frac{T_0(z)}{T_\infty} \right)^{\overline{W}g/\R \Gamma}
   \label{pstrat2}
\ee


\section{Mass and Species Transport}

The species transport equations are solved using a predictor-corrector scheme. Advection terms are written in flux divergence (conservative) form. In the predictor step, the mass density in cell $ijk$ at time level $n+1$ is estimated based on information at the $n$th level:
\be  \frac{(\rho Z_\alpha)_{ijk}^{*}-(\rho Z_\alpha)_{ijk}^n}{\dt}
  + \nabla\!\cdot(\overline{\rho Z_\alpha}^{\rm FL} \mathbf{u})_{ijk}^n
  = \nabla\!\cdot (\rho D_\alpha \nabla Z_\alpha)_{ijk}^n + \left( \dot{m}'''_\alpha + \dot{m}'''_{\rm b,\alpha} \right)_{ijk}^n
\ee
The quantity $\overline{\rho Z_\alpha}^{\rm FL}$ indicates a \emph{flux limiter} applied to the cell face value, as discussed below in Section \ref{sec_flux_limiters}. The mass source terms due to chemistry, evaporation, or pyrolysis are computed at the end of the previous time step and used in both the predictor and corrector steps. The mean chemical source term, $\dm_{\alpha}'''$, is discussed in Chapter~\ref{combustionsection}.  The bulk subgrid source term, $\dm_{\rm b,\alpha}'''$, is discussed in Chapters~\ref{chapter:solid_phase} and \ref{chapter:lagrangian_particles} on solid phase pyrolysis and Lagrangian particles, respectively.

In DNS mode, the molecular diffusivity is based on mixture-averaged binary Fickian diffusion.  In LES mode (default) the diffusivity is taken from the molecular and turbulent viscosities divided by the turbulent Schmidt number.  That is, to save cost we approximate the molecular plus turbulent diffusivity by $(\mu + \mu_t)/\mbox{Sc}_t$.   The turbulent Schmidt number is constant with default value $\mbox{Sc}_t = 0.5$.  The model for the turbulent viscosity $\mu_t$ is discussed in Section \ref{section:turbulent_viscosity}.  Optionally, by setting \textct{SIMULATION\_MODE='LES'} on \textct{MISC}, the molecular and turbulent transport coefficients are treated separately, $\rho D_\alpha + \mu_t/\mbox{Sc}_t$ (at added cost).  The same applies for the thermal diffusivity.

The corrector step is as follows:
\be \frac{(\rho Z_\alpha)_{ijk}^{n+1}-\ha\left[(\rho Z_\alpha)_{ijk}^n
    +(\rho Z_\alpha)_{ijk}^{*}\right]} {\ha \dt}
    + \nabla\!\cdot(\overline{\rho Z_\alpha}^{\rm FL} \mathbf{u})_{ijk}^*
    = \nabla\!\cdot (\rho D_\alpha \nabla Z_\alpha)_{ijk}^{*} + \left( \dot{m}'''_\alpha + \dot{m}'''_{\rm b,\alpha} \right)_{ijk}^n
\ee


\subsection{Flux Limiters}
\label{sec_flux_limiters}

A \emph{flux limiter} is an interpolation scheme for defining mass fluxes at cell faces. Simple linear interpolation of the cell-centered scalar variables to the cell face would result in a central difference scheme.  Such purely centered schemes are known to generate intolerable levels of dispersion error (spurious wiggles) leading to unphysical results such as negative densities or mass fractions outside the range of [0,1].  To address this issue, FDS relies on two schemes: a \emph{flux limiter} (discussed below) that handles the bulk of the problem, and a \emph{flux correction} (see Appendix~\ref{app_boundedness}) that adds the minimum amount of numerical diffusion to maintain boundedness.

For uniform flow velocity, a fundamental property of the exact solution to the equations governing scalar transport is that the total variation of the scalar field (the sum of the absolute values of the scalar differences between neighboring cells) is either preserved or diminished (never increased).  In other words, no new extrema are created.  Numerical schemes which preserve this property are referred to as total variation diminishing (TVD) schemes.  The practical importance of using a TVD scheme for fire modelling is that such a scheme is able to accurately track coherent vortex structure in turbulent flames and does not develop spurious reaction zones.

FDS employs two second-order TVD schemes as options for scalar transport: Superbee and CHARM.  Superbee \cite{Roe:1986} is recommended for LES because it more accurately preserves the scalar variance for coarse grid solutions that are not expected to be smooth.  Due to the gradient steepening applied in Superbee, however, the convergence degrades at small grid spacing for smooth solutions (the method will revert to a stair-step pattern instead of the exact solution).  CHARM \cite{Zhou:1995}, though slightly more dissipative than Superbee, is convergent, and is therefore the better choice for DNS calculations where the flame front is well resolved.

To illustrate how flux limiters are applied to the scalar transport equations, below we discretize the advection terms in Eq.~(\ref{species}) in one dimension:
\be  \frac{(\rho Z)_{i}^* - (\rho Z)_{i}^n}{\dt}
    + \frac{\overline{\rho Z}^{\rm FL}_{i+\frac{1}{2}} u_{i+\frac{1}{2}} - \overline{\rho Z}^{\rm FL}_{i-\frac{1}{2}} u_{i-\frac{1}{2}}}{\dx} = ...
\ee
Note that the $\pm\ha$ suffixes indicates a face value for a particular cell $i$. A flux-limited scalar value (density in this case) premultiplies the staggered, face-centered velocity to form the scalar advective flux.  Recall that these velocity values are primitive variables in the calculation---they are \emph{not} interpolated.
Consider face $i+\frac{1}{2}$ between cells $i$ and $i+1$ and let $\phi$ denote a general scalar variable, like $\rho Z_\alpha$.  The local (loc) and upstream (up) data variations are
\begin{align}
\delta \phi_{\rm loc} &= \phi_{i+1}-\phi_i \\
\delta \phi_{\rm up}  &= \left\{ \begin{array}{ll} \phi_i-\phi_{i-1} & \mbox{if} \quad u_i>0 \\ \phi_{i+2}-\phi_{i+1} & \mbox{if} \quad u_i<0 \end{array} \right.
\end{align}
The limiter function $B(r)$ depends on the upstream-to-local data ratio, $r=\delta \phi_{\rm up}/\delta \phi_{\rm loc}$. In FDS, options for the limiter function include \cite{Toro}:
\begin{table}[H]
\begin{center}
\begin{tabular}{lc}
Flux Limiter                           & $B(r)$                         \\
\hline
Central Difference                     & 1                              \\
Godunov                                & 0                              \\
MINMOD                                 & $\max(0,\min(1,r))$            \\
Superbee \cite{Roe:1986} (LES default) & $\max(0,\min(2r,1),\min(r,2))$ \\
CHARM \cite{Zhou:1995} (DNS default)   & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
MP5 \cite{Suresh:1997}                 & see below
\end{tabular}
\end{center}
\end{table}
\noindent For the Central Difference, Godunov, MINMOD, and Superbee limiters, the scalar face value is found from
\begin{equation}
\label{eqn_flux_limiter}
\overline{\phi}^{\rm FL}_{i+1/2} = \left\{ \begin{array}{lcll} \phi_i &+& B(r) \,\frac{1}{2} \,\delta \phi_{\rm loc} & \mbox{if} \quad u_i>0 \vspace{0.2 cm}\\
\phi_{i+1} &-& B(r) \,\frac{1}{2} \,\delta \phi_{\rm loc} & \mbox{if} \quad u_i<0 \end{array} \right.
\end{equation}
For CHARM, the face value is given by~\cite{Kempf:2003}
\begin{equation}
\label{eqn_charm_limiter}
\overline{\phi}^{\rm FL}_{i+1/2} = \left\{ \begin{array}{lcll} \phi_i &+& B(r) \,\frac{1}{2} \,\delta \phi_{\rm up} & \mbox{if} \quad u_i>0 \vspace{0.2 cm}\\
\phi_{i+1} &-& B(r) \,\frac{1}{2} \,\delta \phi_{\rm up} & \mbox{if} \quad u_i<0 \end{array} \right.
\end{equation}
The MP5 scheme of Suresh and Huynh \cite{Suresh:1997} is based on the keen observation that three points cannot distinguish between extrema and discontinuities.  The functional form of the limiter is not as simple as the three-point schemes described above, so we refer the reader to the original paper or the FDS source code for details.  But the basic idea behind the method is to use a five-point stencil, three upwind and two downwind, to reconstruct the cell face value, considering both accuracy and monotonicity-preserving constraints.  An additional benefit of the MP5 scheme is that it was designed specifically with strong stability-preserving (SSP) Runge-Kutta time discretizations in mind.  The predictor-corrector scheme used by FDS is similar to the second-order SSP scheme described in \cite{Gottlieb:2001}.


\subsubsection{Notes on Implementation}

In practice, we set $r=0$ initially and only compute $r$ if the denominator is not zero.  Note that for $\delta \phi_{loc}=0$, it does not matter which limiter is used: all the limiters yield the same scalar face value.  For CHARM, we set both $r=0$ and $B=0$ initially and only compute $B$ if $r>0$ (this requires data variations to have the same sign). Otherwise, CHARM reduces to Godunov's scheme.

The Central Difference, Godunov, and MINMOD limiters are included for completeness, debugging, and educational purposes.  These schemes have little utility for typical FDS applications.

\subsection{Time Splitting for Mass Source Terms}
\label{sec_time_splitting}

Following the corrector step of the transport scheme, source terms are computed for the next time step.  The source terms are typically related to particle evaporation or combustion, and these processes are computed at the end of the time step. In the case of combustion, the total mass of a grid cell is not changed; rather the species mass fractions change. The mean chemical source term, $\dm_{\alpha}'''$, is discussed in Chapter~\ref{combustionsection}.  The bulk subgrid source term, $\dm_{\rm b,\alpha}'''$, is discussed in Chapters~\ref{chapter:solid_phase} and \ref{chapter:lagrangian_particles} on solid phase pyrolysis and Lagrangian particles, respectively.


\subsection{Boundary Conditions for Temperature, Species Mass Fraction, and Density}
\label{section:TZD_bc}

The gas temperature, species mass fractions, and density are computed at the center of each grid cell. At an exterior boundary, or at
the boundary of an interior obstruction, these values must be computed at the face of the cell that falls at the boundary interface. In general, the temperature at the boundary, $T_{\rm w}$, is computed first, followed by species mass fractions, $Z_{\rm \alpha,w}$, followed by density, $\rho_{\rm w}$. The density is typically determined from the equation of state:
\be  \rho_{\rm w} = \frac{\overline{p}_m}{ \R \, T_{\rm w} \, \sum_\alpha (Z_{\rm \alpha,w}/W_\alpha) }  \ee
Here, $\overline{p}_m$ denotes the background pressure of the gas phase region.

When necessary, the boundary value is linearly extrapolated one half of a grid cell into the ``ghost'' cell for use by the gas phase solver. In the sections below, the value at the center of the gas phase cell adjacent to the boundary is denoted with the subscript ``g'' (for ``gas phase'', \emph{not} ``ghost''), and the value at the boundary by ``w'' (for ``wall'').

\subsubsection{Solid Boundaries}

At a solid boundary, the surface temperature, $T_{\rm w}$, is either specified or computed as described in Chapter~\ref{chapter:solid_phase}. For an LES calculation, the convective heat flux at the surface is determined via an empirical heat transfer coefficient, $h$, and the convective heat flux at the boundary is written:
\be
   k \frac{T_{\rm g} - T_{\rm w}}{\dn/2} = h \; (T_{\rm g}-T_{\rm w})  \label{ebal}
\ee
where $\dn/2$ is the distance between the surface and the center of the adjacent gas phase cell. The convective heat transfer coefficient, $h$, is described in Section~\ref{conflux}. For a DNS calculation, the convective heat transfer is determined directly from the computed or specified surface temperature.

There is no transfer of mass at a solid boundary; thus, the boundary value for the species mixture $\alpha$ is simply
\be Z_{\rm \alpha,w} = Z_{\rm \alpha,g} \ee

\subsubsection{Open Boundaries}

The term ``open'' denotes a non-solid exterior boundary of the computational domain. Gases are allowed to flow freely in and out. At these boundaries, the temperature and species mass fractions take on their respective exterior values if the flow is incoming, and take on their respective values in the grid cell adjacent to the boundary if the flow is outgoing. This is a simple upwind boundary condition.



\subsubsection{Specified Mass Flux}

Here, the mass flux of species $\alpha$, $\dot{m}_\alpha''$, is specified or computed as part of the overall solid phase calculation. To determine the mass fraction of species mixture $\alpha$ at the boundary, $Z_{\alpha,f}$, the following equations must be solved iteratively
\begin{gather}
\label{eqn_total_mass_flux} \sum_\alpha \dot{m}_\alpha'' = \rho_{\rm w} u_n \\
\label{eqn_spec_mass_flux}  \dot{m}_\alpha'' = u_n \rho_{\rm w} Z_{\rm \alpha,w} - (\rho D_\alpha)_{\rm w} \, \frac{Z_{\rm \alpha,g}-Z_{\rm \alpha,w}}{\dn/2}
\end{gather}
where $u_n$ is the normal component of velocity at the wall pointing into the flow domain and $\dn/2$ is the distance between the center of the gas cell and the wall. Together with the equation of state, Eqs.~(\ref{eqn_total_mass_flux}) and (\ref{eqn_spec_mass_flux}) are solved iteratively for the unknowns $\rho_{\rm w}$, $u_n$, and $Z_{\rm \alpha,w}$.  The surface temperature used in the EOS depends on the thermal boundary condition.


\subsubsection{Mesh Interface Boundaries}

In simulations involving more than one numerical mesh, information has to be passed between meshes, even when
the meshes are being processed by separate computers. If two meshes abut each other, and the mesh cells are aligned and the same size, then
one mesh simply uses the density and species mass fractions of the adjacent mesh as the ``ghost'' cell values. However, in cases where the
mesh cells are not the same size, the exchange of information must be done more carefully. Consider a case where two meshes meet:

\begin{figure}[h!]
\begin{picture}(200,110)(0,-10)
\setlength{\unitlength}{0.02in}
\put(120,10){\framebox(20,20){ }}
\put(120,30){\framebox(20,20){ }}
\put(140,10){\framebox(40,40){ }}
\put(100,30){\makebox(0,0){Mesh 1}}
\put(200,30){\makebox(0,0){Mesh 2}}
\thicklines
\put(140,0){\line(0,1){60}}
\end{picture}
\caption{Mesh interface boundary with 2:1 refinement.}
\label{fig:meshinterface}
\end{figure}

\noindent
We want the total and species mass fluxes between meshes to be the same. Let the density in cell $(1,j',k')$ of Mesh 2 be denoted $\rho_{1,j'k'}^{(2)}$. Assume that this cell abuts two cells in Mesh 1. The densities in the two abutting cells of Mesh 1 are denoted $\rho_{I,jk}^{(1)}$. Note that $j$ and $k$ are not the same as $j'$ and $k'$. $I$ is the number of cells in the $x$ direction of Mesh 1. The ghost cell quantities in Mesh 1 have an $i$ index of $I+1$. The ghost cell quantities in Mesh 2 have an $i$ index of 0.
We want to assert mass conservation at the mesh interface:
\be
   \sum_{j,k} u_{I,jk}^{(1)} \; \rho_{{\rm w},jk}^{(1)} \; \dy^{(1)} \, \dz^{(1)}  =
              u_{0,j'k'}^{(2)} \; \rho_{{\rm w},j'k'}^{(2)} \; \dy^{(2)} \, \dz^{(2)}  \label{rhou}
\ee
To enforce this condition, we obtain $\rho_{{\rm w},jk}^{(1)}$ on Mesh 1 and $\rho_{{\rm w},j'k'}^{(2)}$ on Mesh 2 from a flux limiter (see Section \ref{sec_flux_limiters}) once data has been exchanged between meshes.  Since only one layer of ghost cells is exchanged, the second upwind data value is linearly extrapolated to maintain second-order accuracy of the interpolated face density.


\subsubsection{Special Topic: Vertical Mesh Coarsening in Atmospheric Flows}

In atmospheric flows, it may be beneficial to coarsen the mesh in the vertical direction.  Matching the heat flux at the mesh interface is not trivial due to the background pressure stratification and the relative coarse grid resolution usually applied to atmospheric flow calculations.  FDS does not explicitly match the flux.  Instead, it relies on the temperature gradient constructed from the interior ({\ct KKG}) and ghost cell ({\ct KK}) values to be equivalent (additionally, the interpolation of the transport coefficient at the mesh boundary must match, but this is less of a problem than the temperature gradient).  To further complicate matters, FDS always extracts the gas phase temperature value from the ideal gas law, given the cell species composition, density, and background pressure.

\begin{equation}
\label{eq:EOScode}
T(z) = \frac{P(z) \overline{W}(z)}{\rho(z) \, R} = \frac{\mathtt{PBAR}(z)}{\mathtt{RHO}(z) \, \mathtt{RSUM}(z)}
\end{equation}

At {\ct INTERPOATED} boundaries with mesh coarsening (or refinement) the positions of the ghost cell values do not match positions of the interior gas phase cell values from the neighboring mesh.  This is depicted for our problem in Fig.~\ref{fig:vertmeshcoarsening}.  Let $\delta z^{(1)}_{kkg}$ denote the cell height for the interior gas phase cell on Mesh 1, for example.  To match the temperature gradient at the mesh interface we require

\begin{equation}
\label{eq:TMPatminterp}
\frac{T^{(1)}_{kk} - T^{(1)}_{kkg}}{\mbox{$\frac{1}{2}$} \delta z^{(1)}_{kk} + \mbox{$\frac{1}{2}$}\delta z^{(1)}_{kkg}} = \frac{T^{(2)}_{kkg} - T^{(1)}_{kkg}}{\mbox{$\frac{1}{2}$} \delta z^{(2)}_{kkg} + \mbox{$\frac{1}{2}$}\delta z^{(1)}_{kkg}}
\end{equation}

Now, suppose we are filling ghost cell values for Mesh 1.  Our goal is to specify the ghost cell value of the density, $\rho^{(1)}_{kk}$, such that Eq.~(\ref{eq:TMPatminterp}) holds.  Plugging Eq.~(\ref{eq:EOScode}) into Eq.~(\ref{eq:TMPatminterp}), and defining
\begin{equation}
\label{eq:ddofactor}
\mathtt{DDO} = \frac{\delta z^{(1)}_{kk} + \delta z^{(1)}_{kkg}}{\delta z^{(2)}_{kkg} + \delta z^{(1)}_{kkg}}
\end{equation}
we get
\begin{equation}
\label{eq:atmrhoghost}
\rho^{(1)}_{kk} = \left. \left( \frac{P^{(1)}_{kk}}{\mathtt{RSUM}^{(1)}_{kk}} \right) \middle/ \left\{ \left(\frac{P^{(1)}_{kkg}}{\rho^{(1)}_{kkg}\,\mathtt{RSUM}^{(1)}_{kkg}}\right) + \mathtt{DDO}\left[ \left(\frac{P^{(2)}_{kkg}}{\rho^{(2)}_{kkg}\,\mathtt{RSUM}^{(2)}_{kkg}}\right) - \left(\frac{P^{(1)}_{kkg}}{\rho^{(1)}_{kkg}\,\mathtt{RSUM}^{(1)}_{kkg}}\right) \right] \right\} \right.
\end{equation}
A similar expression is obtained when viewing the interpolation from Mesh 2.

Note that this interpolation may be thought of as a central difference for the temperature gradient.  The corresponding density variation is nonlinear.  The mass fluxes computed by the flux limiter scheme described above still match at the refined mesh interface and therefore do not require adjustment.

This special treatment of the ghost cell density may be accessed in the source code by searching on the logical flag {\ct ATMOSPHERIC\_INTERPOLATION}.  By default, the procedure is invoked at vertical mesh refinement boundaries only if {\ct STRATIFICATION=.TRUE.} and if the vertical cell spacing {\ct DZ} on any mesh exceeds 2 m.  For smaller grid spacing the spurious temperature effect is not noticeable and it is deemed more appropriate to keep the temperature variation consistent with the flux limited interface density.  Setting the developer logical flag {\ct USE\_ATMOSPHERIC\_INTERPOLATION} on the {\ct WIND} line overrides the default behavior.

\begin{figure}[h!]

\begin{picture}(200,150)(0,-75)
\setlength{\unitlength}{0.02in}

\linethickness{0.25mm}
\put(140,0){\framebox(40,40){ }}
\put(140,-40){\framebox(40,40){ }}

\linethickness{0.05mm}
\put(140,-20){\framebox(20,20){ }}
\put(160,-20){\framebox(20,20){ }}
\put(140,-40){\framebox(20,20){ }}
\put(160,-40){\framebox(20,20){ }}
\put(140,0){\framebox(20,20){ }}
\put(160,0){\framebox(20,20){ }}

\put(150,10){\circle{2}}
\put(150,-10){\circle*{2}}
\put(170,10){\circle{2}}
\put(170,-10){\circle*{2}}
\put(160,20){\circle*{3}}
\put(160,-20){\circle{3}}

\put(100,-20){\makebox(0,0){Mesh 1}}
\put(100, 20){\makebox(0,0){Mesh 2}}

\put(125,0){\line(1,0){11}}
\put(100,0){\makebox(0,0){Mesh Interface}}

\put(174,-10){\line(1,0){21}}
\put(225,-10){\makebox(0,0){\ct M(1)\%TMP(KKG)}}
\put(253,-10){\vector(1,0){10}}
\put(272,-10){\makebox(0,0){$T^{(1)}_{kkg}$}}
\put(174,10){\line(1,0){21}}
\put(223,10){\makebox(0,0){\ct M(1)\%TMP(KK)}}
\put(253,10){\vector(1,0){10}}
\put(272,10){\makebox(0,0){$T^{(1)}_{kk}$}}

\put(170,-30){\line(1,0){17}}
\put(216,-30){\makebox(0,0){\ct M(2)\%TMP(KK)}}
\put(242,-30){\vector(1,0){10}}
\put(261,-30){\makebox(0,0){$T^{(2)}_{kk}$}}

\put(170,30){\line(1,0){17}}
\put(218,30){\makebox(0,0){\ct M(2)\%TMP(KKG)}}
\put(246,30){\vector(1,0){10}}
\put(265,30){\makebox(0,0){$T^{(2)}_{kkg}$}}

\put(170,-30){\line(-1,1){8}}
\put(170,30){\line(-1,-1){8}}
\end{picture}

\caption[Verticel mesh coarening]{Vertical mesh coarsening.  Mesh 1 is the fine mesh.  Mesh 2 is the coarse mesh.  In this example, the coarsening ratio is 1:2.  Solid dots represent interior gas phase unknown values.  Open circles represent ghost cell values.}
\label{fig:vertmeshcoarsening}
\end{figure}

\section{The Velocity Divergence}

Because of the low Mach number assumption, the velocity divergence (the rate of volumetric expansion) plays an important role in the overall solution scheme.  In the FDS algorithm, the divergence is a surrogate for the energy equation.  The divergence is factored out of the conservative form of the sensible enthalpy equation (see Appendix~\ref{app_divergence}) and when the divergence constraint is satisfied (enforced by the momentum update and solution of the Poisson equation for pressure) the conservative form of the sensible enthalpy equation is satisfied by construction.

For the $m$th zone, with background pressure $\bp_m$, the divergence may be written as
\begin{equation}
\label{eqn_divfromeos}
\nabla\!\cdot \bu = {D} - {P}\; \dod{\bp_m}{t}
\end{equation}
where
\begin{equation}
\label{eqn_fdsP1}
P = \frac{1}{\overline{p}_m} - \frac{1}{\rho c_p T}
\end{equation}
and
\begin{align}
\label{eqn_fdsD1}
D &= \frac{1}{\rho c_p T}\left[ \dot{q}^\ppp + \dot{q}_{\rm b}^\ppp - \Div \dot{\mathbf{q}}^\pp - \mathbf{u} \cdot\nabla (\rho h_{\rm s}) + w \rho_0 g_z\right] \notag\\[.1in]
&+ \frac{1}{\rho} \sum_\alpha \left(\frac{\overline{W}}{W_\alpha} - \frac{h_{\rm s,\alpha}}{c_p T} \right) \bigg[ \Div (\rho D_\alpha \nabla Y_\alpha) - \mathbf{u} \cdot \nabla (\rho Y_\alpha) +\dot{m}_\alpha^\tripleprime  \bigg] \notag \\[.1in]
&+ \frac{1}{\rho} \sum_\alpha \left(\frac{\overline{W}}{W_\alpha} - \frac{\int_{T_{\rm b}}^T c_{p,\alpha}(T') \, {\rm d}T'}{c_p T} \right) \, \dot{m}_{\rm b,\alpha}^\tripleprime
\end{align}

\subsection{Mass and Energy Source Terms}
\label{div_source_terms}

The volumetric source terms in the divergence expression require extended discussion.  The heat release rate per unit volume, $\dq'''$, and the mass generation rate of species $\alpha$ per unit volume, $\dot{m}_\alpha'''$, are detailed in Chapter~\ref{chapter:combustion}, Combustion. The flux term, $\dot{\mathbf{q}}^\pp$, is defined in Eq.~(\ref{bqdot_def}). The radiative source term, $\dq_{\rm r}'''$, that is included in $\nabla \cdot \dot{\mathbf{q}}^\pp$ is discussed in Chapter~\ref{chapter:radiation}, Thermal Radiation.  The bulk heat source from Lagrangian particles, $\dq_{\rm b}'''$, which accounts for convective heat transfer and radiative absorption, is discussed in Chapter~\ref{chapter:lagrangian_particles}, Lagrangian Particles.  The bulk mass source from Lagrangian particles, $\dot{m}_{\rm b,\alpha}'''$, is also found in Chapter~\ref{chapter:lagrangian_particles}.

The source terms are computed in the corrector stage of the time step, following the update of the density and species mass fractions. The terms in Eq.~(\ref{eqn_fdsD1}) involving $\dq_{\rm b}'''$, $\dot{m}_\alpha'''$, and $\dot{m}_{\rm b,\alpha}'''$ are stored in an array called {\ct D\_SOURCE} and applied in the construction of the divergence expression required for the corrected update of the velocity.

\subsection{Diffusion Terms}
\label{div_discret}

The thermal and material diffusion terms of Eq.~(\ref{eqn_fdsD1}) are pure second-order central differences. For example, the thermal
conduction term is differenced as follows:
\begin{align}
(\nabla\!\cdot k \nabla T)_{ijk}
            &&=&&& \frac{1}{\dx} \Bigg[ && k_{i+\ha,jk} && \frac{T_{i+1,jk}-T_{ijk}}{\dx}  &&-&& k_{i-\ha,jk}  && \frac{T_{ijk}-T_{i-1,jk}}{\dx}  &&\Bigg]  &+ \notag \\
            && &&& \frac{1}{\dy} \Bigg[ && k_{i,j+\ha,k}&& \frac{T_{i,j+1,k}-T_{ijk}}{\dy} &&-&& k_{i,j-\ha,k} && \frac{T_{ijk}-T_{i,j-1,k}}{\dy} &&\Bigg]  &+ \notag \\
            && &&& \frac{1}{\dz} \Bigg[ && k_{ij,k+\ha} && \frac{T_{ij,k+1}-T_{ijk}}{\dz}  &&-&& k_{ij,k-\ha}  && \frac{T_{ijk}-T_{ij,k-1}}{\dz}  &&\Bigg]  &
\end{align}
The thermal conductivity at the cell interface, denoted by the $\ha$ cell index, is the average of its values in the two adjacent cells.

\subsection{Corrections for Numerical Mixing}

The differencing of the convection terms, $\mathbf{u} \cdot\nabla (\rho h_s)$ and $\mathbf{u} \cdot \nabla (\rho Y_\alpha)$, is complex.  If not handled carefully, subtle issues related to numerical diffusion in the scalar transport schemes can cause significant conservation errors in the implied energy equation.  The proper discretization of these terms is discussed in Appendix~\ref{app_divergence}.

\subsection{Computing the Temperature}

The mean cell gas temperature, $T$, is derived from the density and species mass fractions via the equation of state:
\be T_{ijk} = \frac{\bp_m}{\rho_{ijk} \R\, \sum_{\alpha=0}^{N_s} (Z_{\alpha,ijk}/W_\alpha)}\ee

\subsection{Sensible Enthalpy}

The sensible enthalpy of the gas is a mass-weighted average of the enthalpies of the lumped species (denoted by $\alpha$), which are in turn a mass-weighted average of the enthalpies of the individual gas species (denoted by $n$):
\be
  h_{\rm s} = \sum_\alpha Z_\alpha \, h_{\rm s,\alpha} \quad;\quad  h_{\rm s,\alpha}=\sum_n Y_n \, h_{{\rm s},n}  \quad; \quad h_{{\rm s},n}(T)=\int_{T_0}^T c_{p,n}(T') \, \mbox{d}T'
\ee
The values of $h_{{\rm s},n}$ and $c_{p,n}$ for the individual gas species are obtained by table lookup from the NIST-JANAF tables~\cite{NIST_JANAF}. The values are taken to the nearest degree Kelvin.

\subsection{Computing the Background Pressure Rise}

To describe how the background pressure of the $m$th pressure zone, $\bp_m$, is updated in time, consider the expression for the
divergence written in compact notation:
\begin{equation}
\label{eqn_divfromeos2}
\nabla\!\cdot \bu = D - P\; \dod{\bp_m}{t}
\end{equation}
The terms $D$ and $P$ are defined by Eqs.~(\ref{eqn_fdsD1}) and (\ref{eqn_fdsP1}), respectively. The subscript $m$ refers to the
number of the {\em pressure zone}; that is, a volume within the computational domain that is allowed to have its own background pressure rise. A closed room
within a building, for example, is a pressure zone.
The time derivative of the background pressure of the $m$th
pressure zone is found by integrating Eq.~(\ref{eqn_divfromeos2}) over the zone volume (denoted by $\Omega_m$):
\begin{equation}
\dod{\bp_m}{t} = \left( \int_{\Omega_m} D \,\d V - \int_{\partial \Omega_m} \bu \cdot \d \bS \right) \Big/ \int_{\Omega_m} P \,\d V  \label{concon2}
\end{equation}
Equation~(\ref{concon2}) is essentially a consistency condition, ensuring that blowing air or starting a fire within a sealed
compartment leads to an appropriate decrease in the divergence within the volume.

\subsection{Combining Pressure Zones}

In the event that a barrier separating two pressure zones should rupture, Eq.~(\ref{concon2}) is modified so that the pressure in the
newly connected zones is driven towards an equilibrium pressure:
\be
  \bp_{\rm eq} = \sum_m \left( \bp_m \int_{\Omega_m} P \, \d V  \right)  \Big/  \sum_m \int_{\Omega_m} P \, \d V \approx \frac{ \sum_m V_m }{ \sum_m (V_m/\bp_m) }
\ee
Note that
\be
  \int_{\Omega_m} P \, \d V \approx  \frac{ V_m}{\gamma \, \bp_m }
\ee
where $V_m$ is the volume of zone $m$ and $\gamma$ is the ratio of specific heats.
To drive the pressure within the connected zones towards each other, a volume flow, $\dot{V}_m^*$, is applied to each zone. This flow is intended to move gas
from zones with the higher pressures towards zones with lower pressures. Eq.~(\ref{concon2}) now becomes:
\be
   \dod{\bp_{\rm eq}}{t} - \frac{ \bp_m - \bp_{\rm eq} }{\tau} =
   \left( \int_{\Omega_m} D \, \d V - \int_{\partial \Omega_m} \bu \cdot \d \bS - \dot{V}_m^* \right) \Big/ \int_{\Omega_m} P \, \d V
\ee
This equation is solved for $\dot{V}_m^*$.
The first term on the left is the change in the equilibrium pressure with time:
\be
   \dod{\bp_{\rm eq}}{t} = \left( \sum_m \int_{\Omega_m} D \, \d V - \sum_m \int_{\partial \Omega_m} \bu \cdot \d \bS \right) \Big/ \sum_m \int_{\Omega_m} P \, \d V
\ee
The summation is over all connected zones, and it is essentially the net change in pressure with time for the entire connected region. If there is any opening to the
exterior of the computational domain, this term is set to zero and all connected zone pressures are driven towards ambient.
The second term on the left forces the pressure in the $m$th pressure zone towards the equilibrium.
The constant, $\tau$, is a characteristic time for the pressure to come into equilibrium. Its default value is on the order of 1~s. In reality, room pressures typically
come into equilibrium very rapidly, but air movements associated with rapid changes in pressure can cause numerical instabilities.


{\bf Note:} Because of the low Mach number assumption, FDS should not be used for rapid discharge of pressure vessels.





